#!/usr/bin/env python
import numpy as np
from astropy.io import fits
import sys
import matplotlib.pyplot as plt
from scipy.stats import sigmaclip

bias_f = sys.argv[1]
mask_f = sys.argv[2]
gain_f = sys.argv[3]
bias = fits.getdata(bias_f)
gain = fits.getdata(gain_f)
#mask = fits.getdata(mask_f)
mask = fits.getdata(mask_f,dtype='>i2')
flags = np.zeros(mask.shape)
bias_clipped,l,h = sigmaclip(bias.flatten(),4,4)
bias_mean = np.mean(bias_clipped)
flags[np.where((bias-bias_mean)>15)] = 1
#flags[np.where(mask<32)] = 0
onum = len(np.where(flags == 1)[0])
print(onum,onum/9216/9232)
fits.writeto("bias_outlier.fits",data=flags,overwrite=True)
